ITS/ITSsddanalysis.C

   1 void ITSsddanalysis (Int_t evNumber1=0,Int_t evNumber2=0)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and fill some histograms.
   7 //
   8 //     Root > .L anal.C   //this loads the macro in memory
   9 //     Root > anal();     //by default process first event
  10 //     Root > anal(2);    //process third event
  11 //Begin_Html
  12 /*
  13 <img src="gif/anal.gif">
  14 */
  15 //End_Html
  16 /////////////////////////////////////////////////////////////////////////
  17
  18 // Dynamically link some shared libs
  19
  20   if (gClassTable->GetID("AliRun") < 0) {
  21     gROOT->LoadMacro("loadlibs.C");
  22     loadlibs();
  23   } else {
  24       delete gAlice;
  25       gAlice=0;
  26   }
  27
  28   // Connect the Root Galice file containing Geometry, Kine and Hits
  29   TString *str = new TString("galice.root");
  30   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(str->Data());
  31   if (!file) file = new TFile(str->Data(),"UPDATE");
  32
  33   // Get AliRun object from file or create it if not on file
  34   //   if (!gAlice) {
  35   gAlice = (AliRun*)file->Get("gAlice");
  36   if (gAlice) printf("AliRun object found on file\n");
  37   if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  38   //}
  39
  40   TH2F *local = new TH2F("local","time vs anode difference",101,-1.01,1.01,210,-210.,210.);
  41   TH2F *local1 = new TH2F("local1","time vs anode difference",101,-5.01,5.01,210,-2100.,2100.);
  42   TH2F *dtah = new TH2F("dtah","anode difference vs drift time (hits)",210,-21000.,21000.,101,-10000.01,10000.01);
  43   TH2F *dtap = new TH2F("dtap","anode difference vs drift time (points)",21000,-21000.,210.,101,-10000.01,10000.01);
  44   TH1F *dclutim = new TH1F("dclutim","cluster time difference (dan < 1)",200,-2000.,2000.);
  45   TH1F *dfkctim = new TH1F("dfkctim","cluster time difference (dan > 5)",200,-2000.,2000.);
  46   TH2F *xy3 = new TH2F("xy3","y vs x",100,-200.02,200.02,100,-200.02,200.02);
  47   TH2F *xz3 = new TH2F("xz3","x vs z",100,-200.02,200.02,100,-200.02,200.02);
  48   TH2F *yz3 = new TH2F("yz3","y vs z",100,-200.02,200.02,100,-200.02,200.02);
  49   TH2F *xy4 = new TH2F("xy4","y vs x",100,-200.02,200.02,100,-200.02,200.02);
  50   TH2F *xz4 = new TH2F("xz4","x vs z",100,-200.02,200.02,100,-200.02,200.02);
  51   TH2F *yz4 = new TH2F("yz4","y vs z",100,-200.02,200.02,100,-200.02,200.02);
  52   TH1F *chd = new TH1F("chargediff","Charge difference (Gen-Rec)",100,-1000.,1000.);
  53   TH1F *chr = new TH1F("chargeratio","Charge ratio (Gen/Rec)",100,0.,0.1);
  54   TH2F *chp = new TH2F("chp","Point Charge vs time",28,0.,7000.,300,0.,3000.);
  55   TH2F *chh = new TH2F("chh","Hit Charge vs time",28,0.,7000.,80,0.,40.);
  56   TH2F *dadth = new TH2F("dadth","danode vs dtime (hits)",560,-14000.,14000.,601,-300.5,300.5);
  57   TH2F *dadt = new TH2F("dadt","danode vs dtime (points)",280,-7000.,7000.,601,-300.5,300.5);
  58   TH2F *aa = new TH2F("aa","anode hit vs point",204,-0.5,50.5.,204,-.5,50.5);
  59   TH2F *at = new TH2F("at","anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
  60   TH2F *tt = new TH2F("tt","time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
  61   TH2F *at1 = new TH2F("at1","(1a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
  62   TH2F *tt1 = new TH2F("tt1","(1a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
  63   TH2F *at2 = new TH2F("at2","(2a) anode difference vs drift path (mm)",18,0.,36.,100,-200.,200.);
  64   TH2F *tt2 = new TH2F("tt2","(2a) time coord. difference (um) vs drift path (mm)",18,0.,36.,100,-200.,200.);
  65   TH2F *asigma = new TH2F("asigma","anode sigma vs drift path (mm)",18,0.,36.,300,0.,300.);
  66   TH2F *tsigma = new TH2F("tsigma","tau. vs drift path (mm)",18,0.,36.,200,0.,100.);
  67   TH2F *asigma2 = new TH2F("asigma2","2 anode sigma vs drift path (mm)",18,0.,36.,150,0.,300.);
  68
  69   TH1F *dtrp = new TH1F("dtrp","Double track separation (microns)",40,0.,2000.);
  70   TH1F *dtrpall = new TH1F("dtrpall","Double track separation (mm)",100,0.,100.);
  71   TH1F *dtrh = new TH1F("dtrh","Double track separation (microns)",40,0.,2000.);
  72   TH1F *dtrhall = new TH1F("dtrhall","Double track separation (mm)",100,0.,100.);
  73
  74   TH1F *p = new TH1F("p","Momentum ",100,0.,20.);
  75
  76   TH1F *effh = new TH1F("effh","Hit Multiplicity vs drift path (mm)",18,0.,36.);
  77   TH1F *effp = new TH1F("effp","Point Multiplicity vs drift path (mm)",18,0.,36.);
  78
  79   TH2F *anodes = new TH2F("nanodes","Anode Multiplicity vs drift time",28,0.,7000.,5,0.5,5.5);
  80   TH2F *andtsm = new TH2F("nand_ntsm","Anode Mult vs Time Mult",15,0.5,15.5,5,0.5,5.5);
  81   TH2F *tsampl = new TH2F("nsampls","Sample Multiplicity vs drift time",28,0.,7000.,15,0.5,15.5);
  82   TH2F *ntotal = new TH2F("ntotal","Cluster Multiplicity vs drift time",28,0.,7000.,60,0.5,60.5);
  83   TH2F *clmult = new TH2F("clmult","Anode Multiplicity vs Total Multiplicity",50,0.5,50.5,5,0.5,5.5);
  84   TH2F *amplit1 = new TH2F("amplit1","Point Amplitude vs drift path (mm)",28,0.,7000.,64,0.5,1024.5);
  85   TH2F *amplit = new TH2F("amplit","Point Amplitude vs drift path (mm)",28,0.,7000.,60,0.5,600.5);
  86   TH1F *hitpnt = new TH1F("hitpnt","Hit-Point Multiplicity",21,-10.5,10.5);
  87
  88   TH1F *nmatch = new TH1F("nmatch","Number of matched points",5,-0.5.,4.5);
  89   TH1F *rec_vs_time = new TH1F("rec_vs_time","Point Rec. vs drift path",36,0.,36.);
  90   TH1F *hit_vs_time = new TH1F("hit_vs_time","Hit vs drift path",36,0.,36.);
  91   TH1F *rec_vs_time3 = new TH1F("rec_vs_time3","Point Rec. vs drift path",36,0.,36.);
  92   TH1F *hit_vs_time3 = new TH1F("hit_vs_time3","Hit vs drift path",36,0.,36.);
  93   TH1F *rec_vs_time4 = new TH1F("rec_vs_time4","Point Rec. vs drift path",36,0.,36.);
  94   TH1F *hit_vs_time4 = new TH1F("hit_vs_time4","Hit vs drift path",36,0.,36.);
  95   TH1F *rec_vs_time1 = new TH1F("rec_vs_time1","Point Rec. vs drift path",36,0.,36.);
  96   TH1F *hit_vs_time1 = new TH1F("hit_vs_time1","Hit vs drift path",36,0.,36.);
  97   TH1F *fake_vs_time = new TH1F("fake_vs_time","fake points vs drift path",36,0.,36.);
  98
  99   TH1F *noihist = new TH1F("noisehist","noise",80,10.,30.);
 100
 101   TH1F *occupancy3 = new TH1F("occupancy3","Occupancy vs Detector Number, Layer 3",20,0.5,20.5);
 102   TH1F *occupancy4 = new TH1F("occupancy4","Occupancy vs Detector Number, Layer 4",20,0.5,20.5);
 103
 104   TH2F *pntmap3 = new TH2F("pntmap3","Point map Layer 3",20,0.5,20.5,10,0.5,10.5);
 105   TH2F *hitmap3 = new TH2F("hitmap3","Hit map Layer 3",20,0.5,20.5,10,0.5,10.5);
 106   TH2F *map3 = new TH2F("map3","Hit/Point map Layer 3",20,0.5,20.5,10,0.5,10.5);
 107   TH2F *pntmap4 = new TH2F("pntmap4","Point map Layer 4",30,0.5,30.5,10,0.5,10.5);
 108   TH2F *hitmap4 = new TH2F("hitmap4","Hit map Layer 4",30,0.5,30.5,10,0.5,10.5);
 109   TH2F *map4 = new TH2F("map4","Hit/Point map Layer 4",30,0.5,30.5,10,0.5,10.5);
 110   TH2F *xz = new TH2F("xz","X vs Z",50,-5,5.,50,-5.,5.);
 111   TH2F *and_tim = new TH2F("and_tim","Tim vs Anode",30,-100,356.,30,-8000.,8000.);
 112   TH2F *pand_ptim = new TH2F("pand_ptim","Tim vs Anode",30,-100,356.,30,-8000.,8000.);
 113
 114   //Int_t nanodes = 256;
 115   //TH2F *mappa3hit[14][6][2];
 116   //TH2F *mappa4hit[22][8][2];
 117   //TH2F *mappa3pnt[14][6][2];
 118   //TH2F *mappa4pnt[22][8][2];
 119   /*
 120     for(Int_t i=0;i<22;i++) {
 121     for(Int_t j=0;j<8;j++) {
 122     for(Int_t k=0;k<2;k++) {
 123     TString *hname = new TString("hitmap_");
 124     TString *cname = new TString("pntmap_");
 125     Char_t lad[2];
 126     sprintf(lad,"%d",i+1);
 127     hname->Append(lad);
 128     hname->Append("_");
 129     cname->Append(lad);
 130     cname->Append("_");
 131     Char_t det[2];
 132     sprintf(det,"%d",j+1);
 133     hname->Append(det);
 134     hname->Append("_");
 135     cname->Append(det);
 136     cname->Append("_");
 137     Char_t wng[2];
 138     sprintf(wng,"%d",k+1);
 139     hname->Append(wng);
 140     cname->Append(wng);
 141     //mappa4hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
 142     //mappa4pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
 143     if(i<14 && j<6) {
 144     mappa3hit[i][j][k] = new TH2F(hname->Data(),hname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
 145     mappa3pnt[i][j][k] = new TH2F(cname->Data(),cname->Data(),nanodes,0.5,nanodes+0.5,256,0.5,256.5);
 146     }
 147     }
 148     }
 149     }
 150   */
 151
 152   AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
 153   if (!ITS) { cout << "no ITS" << endl; return; }
 154
 155   Int_t nparticles = gAlice->GetEvent(0);
 156
 157   Int_t cp[8]={0,0,0,0,0,0,0,0};
 158
 159   AliITSDetType *iDetType=ITS->DetType(1);
 160
 161   AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
 162   if (!res1) {
 163     res1=new AliITSresponseSDD();
 164     ITS->SetResponseModel(1,res1);
 165   }
 166   res1->SetZeroSupp("2D");  // 1D
 167   res1->SetParamOptions("same","same");
 168   //res1->SetFilenames(" ","$(ALICE_ROOT)/ITS/base.dat","$(ALICE_ROOT)/ITS/2D.dat ");
 169   res1->SetCompressParam(cp);
 170   res1->SetDriftSpeed(7.3);
 171   Float_t vdrift = res1->DriftSpeed();
 172
 173   AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
 174   AliITSgeom *aliitsgeo = ITS->GetITSgeom();
 175
 176   Int_t cp[8]={0,0,0,0,0,0,0,0};
 177
 178   Int_t dum = 0;
 179   Float_t apitch = seg1->Dpz(dum);
 180   Float_t tstep = seg1->Dpx(dum);
 181   Float_t maxand = seg1->Npz()/2.;
 182   //cout << "anodes: " << maxand << ", tstep: " << tstep << endl;
 183
 184   Float_t n,b;
 185   res1->GetNoiseParam(n,b);
 186   printf("SDD: noise baseline %f %f zs option %s data type %s\n",n,b,res1->ZeroSuppOption(),res1->DataType());
 187   printf("SDD: DriftSpeed %f TopValue %f\n",res1->DriftSpeed(),res1->DynamicRange());
 188   Float_t dif0,dif1;
 189   res1->DiffCoeff(dif0,dif1);
 190   printf("SDD: dif0 %f dif1 %f\n",dif0,dif1);
 191
 192   AliITSsimulationSDD *sim1=new AliITSsimulationSDD(seg1,res1);
 193   ITS->SetSimulationModel(1,sim1);
 194
 195   //
 196   //   Loop over events
 197   //
 198
 199   Int_t Nh=0;
 200   Int_t Nh1=0;
 201   for (int nev=0; nev<= evNumber2; nev++) {
 202     if(nev>0) {
 203       nparticles = gAlice->GetEvent(nev);
 204       gAlice->SetEvent(nev);
 205     }
 206     cout << "nparticles  " <<nparticles<<endl;
 207     if (nev < evNumber1) continue;
 208     if (nparticles <= 0) return;
 209
 210     // Reset Pointers
 211     AliITShit *itsHit;
 212     AliITSRecPoint *itsPnt = 0;
 213     AliITSRawClusterSDD *itsClu = 0;
 214
 215     // Reset Event Counters
 216     Int_t nGoodTotalHits = 0;
 217     Int_t nGoodTotalPnts = 0;
 218
 219     // Get Hit, Cluster & Recpoints Tree Pointers
 220
 221     TTree *TH = gAlice->TreeH();
 222     Int_t nenthit=TH->GetEntries();
 223     printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
 224
 225     ITS->GetTreeC(nev);
 226     TTree *TC=ITS->TreeC();
 227     Int_t nentclu=TC->GetEntries();
 228     printf("Found %d entries in the Hit tree (must be one per module per event!)\n",nentclu);
 229
 230     TTree *TR = gAlice->TreeR();
 231     Int_t nentrec=TR->GetEntries();
 232     printf("Found %d entries in the Rec tree (must be one per module per event!)\n",nentrec);
 233
 234     // Get Pointers to Clusters & Recpoints TClonesArrays
 235
 236     Int_t iDet = 1;  // 1 = SDD
 237
 238     TClonesArray *ITSclu  = ITS->ClustersAddress(iDet);
 239     TClonesArray *ITSrec  = ITS->RecPoints();
 240
 241     // check recpoints
 242
 243     Int_t nbytes = 0;
 244     Int_t totpoints = 0;
 245     Int_t totclust = 0;
 246
 247     // check hits
 248
 249     Int_t nmodules=0;
 250     ITS->InitModules(-1,nmodules);
 251     ITS->FillModules(nev,0,nmodules,"","");
 252
 253     TObjArray *fITSmodules = ITS->GetModules();
 254
 255     Int_t first = aliitsgeo->GetStartDet(iDet);
 256     Int_t last = aliitsgeo->GetLastDet(iDet);
 257     printf("det type %d first, last %d %d \n",iDet,first,last);
 258
 259     for (Int_t mod=0; mod<last-first+1; mod++) {
 260       cout << "Module: " << mod+1 << endl;
 261       TTree *TR = gAlice->TreeR();
 262       Int_t nentrec=TR->GetEntries();
 263       TClonesArray *ITSrec  = ITS->RecPoints();
 264
 265       ITS->ResetClusters();
 266       TC->GetEvent(mod);
 267       ITS->ResetRecPoints();
 268       nbytes += TR->GetEvent(mod);
 269
 270       Int_t nrecp = ITSrec->GetEntries();
 271       totpoints += nrecp;
 272       //if (nrecp) printf("Found %d rec points for module %d\n",nrecp,mod);
 273       if (!nrecp) continue;
 274
 275       Int_t nrecc = ITSclu->GetEntries();
 276       totclust += nrecc;
 277       //if (nrecc) printf("Found %d clusters for module %d\n",nrecc,mod);
 278
 279       Int_t nrecp = ITSrec->GetEntries();
 280       Int_t startSDD = aliitsgeo->GetStartSDD();
 281       Int_t *flagP = new Int_t [nrecp];
 282       memset( flagP, 0, sizeof(Int_t)*nrecp );
 283
 284       //printf("point loop\n");
 285
 286       Int_t nGoodPnts = 0;
 287       for (Int_t pnt=0;pnt<nrecp;pnt++) {
 288         itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 289         if(!itsPnt) continue;
 290         itsClu  = (AliITSRawClusterSDD*)ITSclu->At(pnt);
 291         if(!itsClu) continue;
 292         //itsClu->PrintInfo();
 293         nGoodPnts++;
 294         nGoodTotalPnts++;
 295
 296         Int_t pntlayer;
 297         Int_t pntladder;
 298         Int_t pntdetector;
 299         aliitsgeo->GetModuleId(mod+first,pntlayer,pntladder,pntdetector);
 300         Int_t pntmult = itsClu->Samples();
 301         Int_t pntands = itsClu->Anodes();
 302         Float_t pnttime = itsClu->T();
 303         Float_t pntanod = itsClu->A();
 304         Float_t pntchrg = itsClu->Q();
 305         Float_t pntampl = itsClu->PeakAmpl();
 306         Float_t pntpath = pnttime*vdrift/1000.;
 307
 308         Float_t wy = 0.;
 309         if(itsClu->Anodes() != 0.) {
 310           wy = pntmult/((Float_t) pntands);
 311         }
 312         clmult->Fill((Float_t)pntmult,(Float_t) pntands);
 313         ntotal->Fill(pnttime,(Float_t)pntmult);
 314         tsampl->Fill(pnttime,wy);
 315         amplit->Fill(pnttime,pntampl);
 316         amplit1->Fill(pnttime,pntampl);
 317
 318         //  Detector  occupancy
 319
 320         if(pntlayer == 3) {
 321           occupancy3->Fill((Float_t) pntdetector,(Float_t) pntmult);
 322           //mappa3pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
 323         }
 324         if(pntlayer == 4) {
 325           occupancy4->Fill((Float_t) pntdetector,(Float_t) pntmult);
 326           //mappa4pnt[pntladder-1][pntdetector-1][0]->Fill(pntanod,pnttime);
 327         }
 328
 329         //  Point  Efficiency  vs  time.
 330
 331         effp->Fill(pntpath);
 332         anodes->Fill(pnttime,pntands);
 333         andtsm->Fill(wy,pntands);
 334         chp->Fill(pnttime,pntchrg);
 335       }
 336
 337       //printf("hit loop\n");
 338
 339       Float_t sddLength = seg1->Dx();
 340       Float_t sddWidth = seg1->Dz();
 341       Float_t driftSpeed=res1->DriftSpeed();
 342
 343       Int_t nGoodHits = 0;
 344
 345       AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod+first);
 346       Int_t nhits = Mod->GetNhits();
 347       for (Int_t hit=0;hit<nhits;hit++) {
 348         itsHit   = (AliITShit*)Mod->GetHit(hit);
 349
 350         Float_t avx = 0.;
 351         Float_t avy = 0.;
 352         Float_t avz = 0.;
 353         Int_t ifl = 0;
 354         Float_t DepEnergy = 100000.*itsHit->GetIonization();
 355         AliITShit *itsHit1 = 0;
 356         if(DepEnergy == 0.) {
 357           hit++;
 358           if(hit == nhits) break;
 359           itsHit1 = (AliITShit*) Mod->GetHit(hit);
 360           avx = itsHit1->GetXG();
 361           avy = itsHit1->GetYG();
 362           avz = itsHit1->GetZG();
 363           ifl = 1;
 364         }
 365         avx += itsHit->GetXG();
 366         avy += itsHit->GetYG();
 367         avz += itsHit->GetZG();
 368         if(DepEnergy == 0.) {
 369           avx /= 2.;
 370           avy /= 2.;
 371           avz /= 2.;
 372         }
 373         if(ifl == 0) continue;
 374
 375         Float_t px; Float_t py; Float_t pz;
 376         itsHit->GetMomentumG(px,py,pz);
 377         Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz);
 378         p->Fill(ptot*100);
 379         if(ptot < 0.05) continue;
 380
 381         Int_t Layer = itsHit->GetLayer();
 382         Int_t Ladder = itsHit->GetLadder();
 383         Int_t Det = itsHit->GetDetector();
 384
 385         Float_t And;
 386         Float_t Tim;
 387         Float_t x = itsHit->GetXL();
 388         Float_t z = itsHit->GetZL();
 389         xz->Fill(z,x);
 390         seg1->GetPadTxz(x,z);
 391         And = z;
 392         Tim = x*tstep;
 393         and_tim->Fill(And,Tim);
 394
 395         Float_t And1;
 396         Float_t Tim1;
 397         Float_t x1;
 398         Float_t z1;
 399         if(itsHit1) {
 400           x1 = itsHit1->GetXL();
 401           z1 = itsHit1->GetZL();
 402           xz->Fill(z1,x1);
 403           seg1->GetPadTxz(x1,z1);
 404           And1 = z1;
 405           Tim1 = x1*tstep;
 406           and_tim->Fill(And1,Tim1);
 407         }
 408         Float_t DepEnergy = 100000.*itsHit->GetIonization();
 409         if(DepEnergy == 0.) DepEnergy = 100000.*itsHit1->GetIonization();
 410         if(DepEnergy < 5.) continue;
 411
 412         if(itsHit1) {
 413           Tim += Tim1;
 414           Tim /= 2.;
 415           And += And1;
 416           And /= 2.;
 417         }
 418         if(And < 0. || And > maxand) { cout << "And: " << And << endl; continue; }
 419         Float_t path = TMath::Abs(Tim)*vdrift/1000.;
 420         hit_vs_time->Fill(path);
 421         if(Layer==3) hit_vs_time3->Fill(path);
 422         if(Layer==4) hit_vs_time4->Fill(path);
 423
 424         nGoodHits++;
 425         nGoodTotalHits++;
 426
 427         //if(Layer == 3) mappa3hit[Ladder-1][Det-1][0]->Fill(And,Tim);
 428         //if(Layer == 4) mappa4hit[Ladder-1][Det-1][0]->Fill(And,Tim);
 429
 430         effh->Fill(path);
 431         Float_t ww = DepEnergy;
 432         chh->Fill(TMath::Abs(Tim),ww);
 433         Int_t inmatches = 0;
 434
 435         Float_t diffmin = 100000.;
 436         Int_t pntmin = -1;
 437
 438         //printf("point loop\n");
 439         for (Int_t pnt=0;pnt<nrecp;pnt++) {
 440           itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 441           if(!itsPnt) continue;
 442           itsClu  = (AliITSRawClusterSDD*)ITSclu->At(pnt);
 443           if(!itsClu) continue;
 444
 445           Int_t LayerP;
 446           Int_t LadderP;
 447           Int_t DetP;
 448           aliitsgeo->GetModuleId(mod+first,LayerP,LadderP,DetP);
 449           Int_t LayerH =  itsHit->GetLayer();
 450           Int_t LadderH =  itsHit->GetLadder();
 451           Int_t DetH =  itsHit->GetDetector();
 452           if(LayerH != LayerP) continue;
 453           if(LadderH != LadderP) continue;
 454           if(DetH != DetP) continue;
 455
 456           Float_t Pand = (Float_t) itsClu->A();
 457           if(Pand < 0 || Pand > maxand) { cout << "Pand: " << Pand << endl; continue; }
 458           Float_t Ptim = (Float_t) itsClu->T();
 459           Float_t Pwng = (Float_t) itsClu->W();
 460           if(Pwng == 1) Ptim *= -1.;
 461           pand_ptim->Fill(Pand,Ptim);
 462           Float_t adiff = And-Pand;
 463           Float_t tdiff = Tim-Ptim;
 464           if(And < 0) {
 465             printf("tim %f\n",Tim);
 466             printf("and %f\n",And);
 467           }
 468           if(Pwng == 1) tdiff *=-1.;
 469           local1->Fill(adiff,tdiff);
 470
 471           if(TMath::Abs(adiff) >= 1) continue;
 472           if(TMath::Abs(tdiff) >= 100) continue;
 473
 474           Float_t apdiff = adiff*apitch;
 475           Float_t tpdiff = tdiff*vdrift;
 476
 477           Float_t diff = TMath::Sqrt( apdiff*apdiff+tpdiff*tpdiff );
 478           if( diff < diffmin ){
 479             diffmin = diff;
 480             pntmin = pnt;
 481           }
 482
 483           if(TMath::Abs(adiff) < 1. && TMath::Abs(tdiff) < 100.) {
 484             inmatches++;
 485           }
 486         } // loop (points)
 487
 488         if( pntmin > -1 ) {
 489           if( flagP[pntmin] == 1) continue;
 490           flagP[pntmin] = 1;
 491           itsClu  = (AliITSRawClusterSDD*)ITSclu->At( pntmin );
 492           Float_t Pand = (Float_t) itsClu->A();
 493           Float_t Ptim = (Float_t) itsClu->T();
 494           Float_t Pwng = (Float_t) itsClu->W();
 495           Float_t sigma = itsClu->Asigma();
 496           Float_t tau = itsClu->Tsigma();
 497           Int_t pntands = itsClu->Anodes();
 498           if(Pwng == 1) Ptim *= -1.;
 499           Float_t adiff = And-Pand;
 500           Float_t tdiff = Tim-Ptim;
 501           if(Pwng == 1) tdiff *=-1.;
 502           local->Fill(adiff,tdiff);
 503
 504           Float_t dpath = Ptim*vdrift/1000.;
 505           Float_t dpathh = Tim*vdrift/1000.;
 506           Float_t adpath = TMath::Abs(dpath);
 507           Float_t adpathh = TMath::Abs(dpathh);
 508           Float_t apdiff = adiff*apitch;
 509           Float_t tpdiff = tdiff*vdrift;
 510           aa->Fill(Pand,And);
 511
 512           Int_t pntands = itsClu->Anodes();
 513           if(pntands == 1) {
 514             at1->Fill(adpath,apdiff);
 515             tt1->Fill(adpath,tpdiff);
 516           }
 517           if(pntands == 2) {
 518             at2->Fill(adpath,apdiff);
 519             tt2->Fill(adpath,tpdiff);
 520           }
 521
 522           at->Fill(adpathh,apdiff);
 523           tt->Fill(adpathh,tpdiff);
 524           asigma->Fill(adpathh,sigma);
 525           tsigma->Fill(adpathh,tau);
 526           if(pntands == 2) asigma2->Fill(adpathh,sigma);
 527
 528           Float_t *lP = new Float_t[3];
 529           lP[0] = itsPnt->GetX();
 530           lP[1] = 0.;
 531           lP[2] = itsPnt->GetZ();
 532           Float_t *gP = new Float_t[3];
 533           aliitsgeo->LtoG(LayerH,LadderH,DetH,lP,gP);
 534           Float_t dx = avx - gP[0];
 535           Float_t dy = avy - gP[1];
 536           Float_t dz = avz - gP[2];
 537           delete lP;
 538           delete gP;
 539
 540           Float_t pntchrg = itsClu->Q();
 541           Float_t dq = DepEnergy/0.122 - pntchrg;
 542           Float_t rq = 0;
 543           if(pntchrg != 0) rq = DepEnergy/0.122/((Float_t) pntchrg);
 544           if(LayerH == 3) {
 545             xy3->Fill(dx,dy);
 546             xz3->Fill(dz,dx);
 547             yz3->Fill(dz,dy);
 548           } else if(LayerH == 4) {
 549             xy4->Fill(dx,dy);
 550             xz4->Fill(dz,dx);
 551             yz4->Fill(dz,dy);
 552           }
 553           chd->Fill(dq);
 554           if(rq != 0.) chr->Fill(rq);
 555
 556           rec_vs_time->Fill(adpathh);
 557           if(Layer==3) rec_vs_time3->Fill(adpathh);
 558           if(Layer==4) rec_vs_time4->Fill(adpathh);
 559         }
 560
 561         nmatch->Fill(inmatches);
 562       }  // loop hits
 563
 564       if(nGoodHits != nGoodPnts) {
 565         printf("module: %d",mod+1);
 566         printf(", nGoodHits: %d",nGoodHits);
 567         printf(", nGoodPnts: %d\n",nGoodPnts);
 568       }
 569       Float_t nHP = (Float_t) nGoodHits-nGoodPnts;
 570       hitpnt->Fill(nHP);
 571
 572       Int_t www = 0.;
 573       if(nGoodHits != 0) www = nGoodPnts/((Float_t) nGoodHits);
 574       if(Layer == 3) {
 575         pntmap3->Fill(Ladder,Det,(Float_t) nGoodPnts);
 576         hitmap3->Fill(Ladder,Det,(Float_t) nGoodHits);
 577         map3->Fill(Ladder,Det,www);
 578       }
 579       if(Layer == 4) {
 580         pntmap4->Fill(Ladder,Det,(Float_t) nGoodPnts);
 581         hitmap4->Fill(Ladder,Det,(Float_t) nGoodHits);
 582         map4->Fill(Ladder,Det,www);
 583       }
 584
 585       //printf("double hit loop\n");
 586
 587       Stat_t wh = 1.;
 588       if(nGoodHits > 1)
 589         wh /= (((Float_t) nGoodHits)*((Float_t) nGoodHits)-1)/2.;
 590       else
 591         wh = 0.;
 592
 593       Int_t *flag = new Int_t[nhits];
 594       Int_t nGoodHitsOK = 0;
 595       for (Int_t hit=0;hit<nhits;hit++) {
 596         flag[hit] = 0;
 597         itsHit   = (AliITShit*)Mod->GetHit(hit);
 598         Float_t avx = 0.;
 599         Float_t avy = 0.;
 600         Float_t avz = 0.;
 601         Int_t ifl = 0;
 602         Float_t DepEnergy = 100000.*itsHit->GetIonization();
 603         AliITShit *itsHit1 = 0;
 604         if(DepEnergy == 0.) {
 605           hit++;
 606           flag[hit] = 0;
 607           if(hit == nhits) break;
 608           itsHit1 = (AliITShit*) Mod->GetHit(hit);
 609           avx = itsHit1->GetXG();
 610           avy = itsHit1->GetYG();
 611           avz = itsHit1->GetZG();
 612           ifl = 1;
 613         }
 614         avx += itsHit->GetXG();
 615         avy += itsHit->GetYG();
 616         avz += itsHit->GetZG();
 617         if(DepEnergy == 0.) {
 618           avx /= 2.;
 619           avy /= 2.;
 620           avz /= 2.;
 621         }
 622         if(DepEnergy < 5. && DepEnergy > 0.) continue;
 623         if(ifl == 0) continue;
 624
 625         Float_t px; Float_t py; Float_t pz;
 626         itsHit->GetMomentumG(px,py,pz);
 627         Float_t ptot = TMath::Sqrt(px*px+py*py+pz*pz);
 628         if(ptot < 0.05) continue;
 629
 630         for (Int_t hit1=hit+1;hit1<nhits;hit1++) {
 631           itsHit2   = (AliITShit*)Mod->GetHit(hit1);
 632
 633           Float_t avx2 = 0.;
 634           Float_t avy2 = 0.;
 635           Float_t avz2 = 0.;
 636           Int_t ifl2 = 0;
 637           Float_t DepEnergy2 = 100000.*itsHit2->GetIonization();
 638           AliITShit *itsHit3 = 0;
 639           if(DepEnergy2 == 0.) {
 640             hit1++;
 641             itsHit3 = (AliITShit*) Mod->GetHit(hit1);
 642             avx2 = itsHit3->GetXG();
 643             avy2 = itsHit3->GetYG();
 644             avz2 = itsHit3->GetZG();
 645             ifl2 = 1;
 646           }
 647           avx2 += itsHit2->GetXG();
 648           avy2 += itsHit2->GetYG();
 649           avz2 += itsHit2->GetZG();
 650           if(DepEnergy2 == 0.) {
 651             avx2 /= 2.;
 652             avy2 /= 2.;
 653             avz2 /= 2.;
 654           }
 655           if(DepEnergy2 < 5. && DepEnergy2 > 0.) continue;
 656           if(itsHit->GetLayer() != itsHit2->GetLayer()) continue;
 657           if(itsHit->GetLadder() != itsHit2->GetLadder()) continue;
 658           if(itsHit->GetDetector() != itsHit2->GetDetector()) continue;
 659           if(ifl2 == 0) continue;
 660
 661           Float_t px1; Float_t py1; Float_t pz1;
 662           itsHit2->GetMomentumG(px1,py1,pz1);
 663           Float_t ptot1 = TMath::Sqrt(px1*px1+py1*py1+pz1*pz1);
 664           if(ptot1 < 0.05) continue;
 665
 666           Float_t And;
 667           Float_t Tim;
 668           Float_t x = itsHit->GetXL();
 669           Float_t z = itsHit->GetZL();
 670           seg1->GetPadTxz(x,z);
 671           And = z;
 672           Tim = x*tstep;
 673           if(And < 0 || And > maxand) continue;
 674           Float_t And2;
 675           Float_t Tim2;
 676           Float_t x2 = itsHit2->GetXL();
 677           Float_t z2 = itsHit2->GetZL();
 678           seg1->GetPadTxz(x2,z2);
 679           And2 = z2;
 680           Tim2 = x2*tstep;
 681           if(And2 < 0 || And2 > maxand) continue;
 682           Float_t da = apitch*(And-And2);
 683           Float_t dt = vdrift*(Tim-Tim2);
 684           Float_t danh = And-And2;
 685           Float_t dtmh = Tim-Tim2;
 686           Float_t dist = TMath::Sqrt(da*da+dt*dt);
 687           if(dt < 1000.) {
 688             Float_t wx = dt*clock/(1000.*vdrift);
 689             Float_t wy = da/apitch;
 690             dtah->Fill(wx,wy);
 691           }
 692           if(dist<20.) { cout << "skip hit " << hit1 << endl; flag[hit] = 1; }
 693
 694           if(ifl == 1 && ifl2 == 1) {
 695             if(dist>10)dtrh->Fill(dist,wh);
 696             dtrhall->Fill(dist/1000.,wh);
 697             dadth->Fill(dtmh,danh);
 698           }
 699         }  // end cluster loop
 700
 701         Float_t path = TMath::Abs(Tim)*vdrift/1000.;
 702         if(flag[hit] == 0) { nGoodHitsOK++; hit_vs_time1->Fill(path); }
 703       } // end hit loop
 704       delete [] flag;
 705       printf("nGoodHits: %d",nGoodHits);
 706       printf(", nGoodHitsOK: %d\n",nGoodHitsOK);
 707
 708       //printf("cluster loop\n");
 709
 710       AliITSRawClusterSDD *itsCluFake = 0;
 711       Int_t nGoodPntsOK = 0;
 712       for( int ip=0; ip<nrecp; ip++) {
 713         itsClu  = (AliITSRawClusterSDD*)ITSclu->At(ip);
 714         if(!itsClu) continue;
 715         Float_t Ptim = (Float_t) itsClu->T();
 716         Float_t dpath = Ptim*vdrift/1000.;
 717         Float_t adpath = TMath::Abs(dpath);
 718         if( flagP[ip] == 1) { nGoodPntsOK++; rec_vs_time1->Fill(dpath); }
 719         else {
 720           cout << "ip: " << ip << endl;
 721           itsCluFake  = (AliITSRawClusterSDD*)ITSclu->At( ip );
 722           if(!itsCluFake) continue;
 723           cout << "pointer: " << itsCluFake << endl;
 724           itsCluFake->PrintInfo();
 725           Float_t Ptim = (Float_t) itsCluFake->T();
 726           Float_t dpath = Ptim*vdrift/1000.;
 727           fake_vs_time->Fill(dpath);
 728         }
 729       }
 730
 731       Stat_t wp = 1.;
 732       if(nGoodPntsOK > 1)
 733         wp /= (((Float_t) nGoodPntsOK)*((Float_t) nGoodPntsOK)-1)/2.;
 734       else
 735         wp = 0.;
 736
 737       //printf("double cluster loop\n");
 738       for (Int_t pnt=0;pnt<nrecp;pnt++) {
 739         if( flagP[pnt] == 0) continue;
 740         itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 741         if(!itsPnt) continue;
 742         itsClu  = (AliITSRawClusterSDD*)ITSclu->At(pnt);
 743         if(!itsClu) continue;
 744         AliITSRecPoint *itsPnt1;
 745         AliITSRawClusterSDD *itsClu1;
 746         for (Int_t pnt1=pnt+1;pnt1<nrecp;pnt1++) {
 747           if( flagP[pnt1] == 0) continue;
 748           itsPnt1 = (AliITSRecPoint*)ITSrec->At(pnt1);
 749           if(!itsPnt1) continue;
 750           itsClu1 = (AliITSRawClusterSDD*)ITSclu->At(pnt1);
 751           if(!itsClu1) continue;
 752
 753           Float_t dan = itsClu->A()-itsClu1->A();
 754           Float_t Pwng = (Float_t) itsClu->W();
 755           Float_t dt1 = itsClu->T();
 756           if(Pwng == 1) dt1 *= -1.;
 757           Float_t dt2 = itsClu1->T();
 758           Float_t Pwng = (Float_t) itsClu1->W();
 759           if(Pwng == 1) dt2 *= -1.;
 760           Float_t dtm = dt1-dt2;
 761           dadt->Fill(dtm,dan);
 762           Float_t dap = apitch*(itsClu->A()-itsClu1->A());
 763           Float_t dtp = vdrift*(dt1-dt2);
 764           Float_t distp = TMath::Sqrt(dap*dap+dtp*dtp);
 765           if(TMath::Abs(dan) < 1.) dclutim->Fill(dtp);
 766           if(TMath::Abs(dan) > 10.) dfkctim->Fill(dtp);
 767           if(dtp < 1000.) {
 768             Float_t wx = dtp*clock;
 769             Float_t wy = dap/apitch;
 770             dtap->Fill(wx,wy);
 771           }
 772           dtrp->Fill(distp,wp);
 773           dtrpall->Fill(distp/1000.,wp);
 774         } // end loop points
 775       } // end loop points
 776
 777       delete [] flagP;
 778     } // end loop modules
 779
 780     ITS->ClearModules();
 781     gAlice->CleanDetectors();
 782
 783   } // end loop events
 784
 785   cout << "open output file" << endl;
 786   TFile fhistos("SDD_histos_test.root","RECREATE");
 787   local->Write();
 788   local1->Write();
 789
 790   aa->Write();
 791   at->Write();
 792   tt->Write();
 793   at1->Write();
 794   tt1->Write();
 795   at2->Write();
 796   tt2->Write();
 797   asigma->Write();
 798   tsigma->Write();
 799   asigma2->Write();
 800   dadt->Write();
 801   dadth->Write();
 802   dclutim->Write();
 803   dfkctim->Write();
 804   xy3->Write();
 805   xz3->Write();
 806   yz3->Write();
 807   xy4->Write();
 808   xz4->Write();
 809   yz4->Write();
 810   chd->Write();
 811   chr->Write();
 812   chh->Write();
 813   chp->Write();
 814   dtrp->Write();
 815   dtrpall->Write();
 816   dtah->Write();
 817   dtap->Write();
 818   dtrh->Write();
 819   dtrhall->Write();
 820   effh->Write();
 821   effp->Write();
 822   rec_vs_time->Write();
 823   hit_vs_time->Write();
 824   hit_vs_time1->Write();
 825   rec_vs_time1->Write();
 826   rec_vs_time3->Write();
 827   hit_vs_time3->Write();
 828   rec_vs_time4->Write();
 829   hit_vs_time4->Write();
 830   fake_vs_time->Write();
 831   p->Write();
 832   nmatch->Write();
 833   anodes->Write();
 834   tsampl->Write();
 835   ntotal->Write();
 836   clmult->Write();
 837   andtsm->Write();
 838   amplit->Write();
 839   amplit1->Write();
 840   hitpnt->Write();
 841   noihist->Write();
 842   occupancy3->Write();
 843   occupancy4->Write();
 844   map3->Write();
 845   hitmap3->Write();
 846   pntmap3->Write();
 847   map4->Write();
 848   hitmap4->Write();
 849   pntmap4->Write();
 850   /*
 851     for(Int_t i=0;i<22;i++) {
 852     for(Int_t j=0;j<8;j++) {
 853     for(Int_t k=0;k<2;k++) {
 854     //mappa4hit[i][j][k]->Write();
 855     //mappa4pnt[i][j][k]->Write();
 856     if(i<14 && j<6) {
 857     //mappa3hit[i][j][k]->Write();
 858     //mappa3pnt[i][j][k]->Write();
 859     }
 860     }
 861     }
 862     }
 863   */
 864   xz->Write();
 865   and_tim->Write();
 866   pand_ptim->Write();
 867   fhistos.Close();
 868   file->Close();
 869 }
 870
 871
 872
 873
 874
 875