ITS/dEdXgeant.C

   1
   2 Int_t max_modules=2269;
   3
   4 void dEdXyy (Int_t evNumber1=0,Int_t evNumber2=0,AliITSPid* pid=0,int mtest=0)
   5 {
   6    if (gClassTable->GetID("AliRun") < 0) {
   7       gROOT->LoadMacro("loadlibs.C");
   8       loadlibs();
   9    } else {
  10       delete gAlice;
  11       gAlice=0;
  12    }
  13
  14 // Connect the Root Galice file containing Geometry, Kine and Hits
  15
  16    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  17    printf("file %p\n",file);
  18    if (file) file->Close();
  19    file = new TFile("galice.root","UPDATE");
  20    file->ls();
  21
  22    printf ("I'm after Map \n");
  23
  24 // Get AliRun object from file or create it if not on file
  25
  26    if (!gAlice) {
  27        gAlice = (AliRun*)file->Get("gAlice");
  28        if (gAlice) printf("AliRun object found on file\n");
  29        if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  30    }
  31    printf ("I'm after gAlice \n");
  32
  33  // Create Histogramms -------------------------------------------------
  34
  35      TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,400.);
  36      //TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,400.);
  37      TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,3000.);
  38      //TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,6.);
  39      TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,30.);
  40      //TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,2000.);
  41      TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,15000.);
  42      //TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,6.);
  43      TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,30.);
  44      TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,400.);
  45      //TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,400.);
  46      TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,3000.);
  47      //TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,6.);
  48      TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,30.);
  49      //TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,200.);
  50      TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.);
  51      //TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,6.);
  52      TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,30.);
  53
  54
  55
  56      /*
  57      TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,4000.);
  58      TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,4000.);
  59      TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,40.);
  60      TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,20000.);
  61      TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,40.);
  62      TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,4000.);
  63      TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,4000.);
  64      TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,40.);
  65      TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.);
  66      TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,40.);
  67      */
  68
  69      TH1F *pathSDD = new TH1F("pathSDD","Path length in the SDD",100,0.,1000.);
  70      TH1F *pathSSD = new TH1F("pathSSD","Path length in the SSD",100,0.,1000.);
  71
  72  // -------------------------------------------------------------
  73    AliITS *ITS  = (AliITS*) gAlice->GetModule("ITS");
  74    if (!ITS) return;
  75    AliITSgeom *geom = ITS->GetITSgeom();
  76
  77    // Set the models for cluster finding
  78
  79    // SPD
  80
  81    ITS->MakeTreeC();
  82    Int_t nparticles=gAlice->GetEvent(0);
  83
  84    AliITSDetType *iDetType=ITS->DetType(0);
  85    AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
  86    TClonesArray *dig0  = ITS->DigitsAddress(0);
  87    TClonesArray *recp0  = ITS->ClustersAddress(0);
  88    AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0);
  89    ITS->SetReconstructionModel(0,rec0);
  90
  91    // SDD
  92
  93    AliITSDetType *iDetType=ITS->DetType(1);
  94    AliITSgeom *geom = ITS->GetITSgeom();
  95
  96    AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
  97    if (!seg1) seg1 = new AliITSsegmentationSDD(geom);
  98    AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
  99    if (!res1) res1=new AliITSresponseSDD();
 100
 101    Float_t baseline,noise;
 102    res1->GetNoiseParam(noise,baseline);
 103    Float_t noise_after_el = res1->GetNoiseAfterElectronics();
 104    Float_t thres = baseline;
 105    thres += (4.*noise_after_el);  // TB // (4.*noise_after_el);
 106    printf("thres %f\n",thres);
 107    //res1->Print();
 108
 109    TClonesArray *dig1  = ITS->DigitsAddress(1);
 110    TClonesArray *recp1  = ITS->ClustersAddress(1);
 111    AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1);
 112    rec1->SetCutAmplitude((int)thres);
 113    ITS->SetReconstructionModel(1,rec1);
 114    //rec1->Print();
 115
 116    // SSD
 117
 118    AliITSDetType *iDetType=ITS->DetType(2);
 119    AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
 120    TClonesArray *dig2  = ITS->DigitsAddress(2);
 121    AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
 122    ITS->SetReconstructionModel(2,rec2);
 123
 124 //
 125 //   Loop over events
 126 //
 127    Int_t Nh=0;
 128    Int_t Nh1=0;
 129    for (int nev=0; nev<= evNumber2; nev++) {
 130      Int_t nparticles = 0;
 131      nparticles = gAlice->GetEvent(nev);
 132      cout << "nev         " << nev <<endl;
 133      cout << "nparticles  " << nparticles <<endl;
 134      if (nev < evNumber1) continue;
 135      if (nparticles <= 0) return;
 136
 137      AliITShit *itsHit;
 138      AliITShit *itsHitPrev;
 139      AliITSRecPoint *itsPnt = 0;
 140      AliITSRawClusterSDD *itsClu = 0;
 141
 142      // Get Hit, Cluster & Recpoints Tree Pointers
 143
 144      TTree *TH = gAlice->TreeH();
 145      Int_t nenthit=TH->GetEntries();
 146      printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
 147
 148      ITS->GetTreeC(nev);
 149      TTree *TC=ITS->TreeC();
 150      Int_t nentclu=TC->GetEntries();
 151      printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
 152
 153      TTree *TR = gAlice->TreeR();
 154      Int_t nentrec=TR->GetEntries();
 155      printf("Found %d entries in the RecPoints tree\n",nentrec);
 156
 157      // Get Pointers to Clusters & Recpoints TClonesArrays
 158
 159      TClonesArray *ITSclu  = ITS->ClustersAddress(1);
 160      printf ("ITSclu %p \n",ITSclu);
 161      //            Int_t ncl = ITSclu->GetEntries();
 162      //     cout<<"ncluster ="<<ncl<<endl;
 163      TClonesArray *ITSrec  = ITS->RecPoints();
 164      printf ("ITSrec %p \n",ITSrec);
 165
 166      // check recpoints
 167
 168      //Int_t nbytes = 0;
 169      Int_t totpoints = 0;
 170      Int_t totclust = 0;
 171
 172      // check hits
 173
 174      Int_t nmodules=0;
 175
 176      ITS->InitModules(-1,nmodules);
 177      ITS->FillModules(nev,0,nmodules,"","");
 178
 179      TObjArray *fITSmodules = ITS->GetModules();
 180
 181      Int_t first0 = geom->GetStartDet(0);  // SPD
 182      Int_t last0 = geom->GetLastDet(0);    // SPD
 183      Int_t first1 = geom->GetStartDet(1);  // SDD
 184      Int_t last1 = geom->GetLastDet(1);    // SDD
 185      Int_t first2 = geom->GetStartDet(2);  // SSD
 186      Int_t last2 = geom->GetLastDet(2);    // SSD
 187
 188      //  For the SPD: first0 = 0, last0 = 239     (240 modules);
 189      //  for the SDD: first1 = 240, last1 = 499   (260 modules);
 190      //  for the SSD: first2 = 500, last2 = 2269  (1770 modules).
 191
 192      Int_t mod;
 193      printf("det type %d first0, last0 %d %d \n",0,first0,last0);
 194      printf("det type %d first1, last1 %d %d \n",1,first1,last1);
 195      printf("det type %d first2, last2 %d %d \n",2,first2,last2);
 196        Int_t negtrSDD = 0;
 197        Int_t negtrSSD = 0;
 198
 199     for (mod=0; mod<last2+1; mod++) {
 200       if(mod>max_modules)continue;
 201       if(mod < first1) continue;  // for the SDD/SSD only
 202
 203        AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
 204        Int_t nhits = Mod->GetNhits();
 205        Int_t layer = Mod->GetLayer();
 206        Int_t ladder = Mod->GetLadder();
 207        Int_t detector = Mod->GetDet();
 208        Int_t hit;
 209        Int_t parent;
 210        Int_t dray;
 211        Int_t hitstat;
 212        Int_t partcode;
 213        Int_t hitprim;
 214        Float_t epart;
 215        Float_t pathInSi=300;
 216        Float_t xhit;
 217        Float_t yhit;
 218        Float_t zhit;
 219        Float_t px;
 220        Float_t py;
 221        Float_t pz;
 222        Float_t pmod;
 223        Float_t xhit0 = 1e+7;
 224        Float_t yhit0 = 1e+7;
 225        Float_t zhit0 = 1e+7;
 226
 227        TTree *TR = gAlice->TreeR();
 228        ITS->ResetClusters();
 229        //cout << "CLUSTERS: get" << endl;
 230        TC->GetEvent(mod);
 231        //cout << "RECPOINTS: reset" << endl;
 232        ITS->ResetRecPoints();
 233        //cout << "RECPOINTS: get" << endl;
 234        //TR->GetEvent(mod+1);
 235        TR->GetEvent(mod);
 236
 237        Int_t nrecp = ITSrec->GetEntries();
 238        totpoints += nrecp;
 239        if (!nrecp) continue;
 240
 241        Int_t trackRecp[3];
 242        Int_t itr;
 243        Int_t TrackRecPoint;
 244        Float_t PmodRecP;
 245        Float_t signalRP;
 246
 247 //         cout <<" module,layer,ladder,detector,nrecp,nhits ="<<
 248 //       mod<<","<<layer<<","<<ladder<<","<<detector<<","<<nrecp<<","<<nhits<< endl;
 249        cout<<".";
 250
 251        // ---------- RecPoint signal analysis for the SDD/SSD  ---------
 252
 253        for (Int_t pnt=0;pnt<nrecp;pnt++) {          // recpoint loop
 254
 255          itsPnt  = (AliITSRecPoint*)ITSrec->At(pnt);
 256
 257          Int_t RecPointPrim = 0;
 258          if(!itsPnt) continue;
 259
 260            signalRP = itsPnt->GetQ();
 261            trackRecp[0] = itsPnt->GetLabel(0);
 262            trackRecp[1] = itsPnt->GetLabel(1);
 263            trackRecp[2] = itsPnt->GetLabel(2);
 264
 265 //         cout<<"New Recp: pnt,signal,tr0,tr1,tr2 ="<<pnt
 266 //             <<","<<signalRP<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl;
 267
 268            /*
 269             if(trackRecp[0]<0&&trackRecp[1]<0&&trackRecp[2]<0) {
 270 cout<<"pnt,tr0,1,2 ="<<pnt<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl;
 271  if(mod<first2) negtrSDD += 1;
 272  if(mod>=first2) negtrSSD += 1;
 273             }
 274            */
 275
 276           xhit0 = 1e+7;
 277           yhit0 = 1e+7;
 278           zhit0 = 1e+7;
 279
 280        // Hit loop
 281         for (hit=0;hit<nhits;hit++) {
 282
 283          itsHit   = (AliITShit*)Mod->GetHit(hit);
 284
 285          hitprim = 0;
 286          dray = 0;
 287          Int_t hitRecp = 0;
 288          Int_t trackHit = itsHit->GetTrack();
 289          hitstat = itsHit->GetTrackStatus();
 290          zhit = 10000*itsHit->GetZL();
 291          xhit = 10000*itsHit->GetXL();
 292          yhit = 10000*itsHit->GetYL();
 293          Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
 294
 295          //      cout<<" New hit: hit,track,hitstat,yhit,ehit ="<<hit<<","<<trackHit<<","<<hitstat<<","<<yhit<<","<<dEn<<endl;
 296
 297          // check the primary particle
 298          partcode = itsHit->GetParticle()->GetPdgCode();
 299          parent = itsHit->GetParticle()->GetFirstMother();
 300          if(parent < 0) hitprim = 1; // primery particle
 301
 302
 303
 304          // select the primery hits with the track number equaled to
 305          // the one of the RecPoint
 306          for(itr=0;itr<3;itr++) {
 307            if(trackRecp[itr] == trackHit) hitRecp = 1;
 308            if(trackRecp[itr] == trackHit && hitprim == 1) {
 309              RecPointPrim = 1;
 310              TrackRecPoint = trackRecp[itr];
 311            }
 312            if(trackRecp[itr] == trackHit && hitprim == 0) trackRecp[itr]=-100;
 313          }
 314
 315          if(hitRecp != 1) continue; // hit doesn't correspond to the recpoint
 316
 317          //cout<<"Hit Ok: trackRecp0,1,2,hitRecp,RecPointPrim,TrackRecPoint ="<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<","<<hitRecp<<","<<RecPointPrim<<","<<TrackRecPoint<<endl;
 318
 319          //      if(hitstat == 66 && yhit < -136.) {  // entering hit
 320          if(hitstat == 66 && hitprim == 1) {  // entering hit
 321             xhit0 = xhit;
 322             yhit0 = yhit;
 323             zhit0 = zhit;
 324          }
 325          //              cout<<" xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl;
 326
 327           if(hitstat == 66) continue; // Take only the hits with the not
 328                                       // zero energy
 329
 330           if(xhit0 > 9e+6 || zhit0 > 9e+6) {
 331             // cout<<"default xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl;
 332             continue;
 333           }
 334
 335          // check the particle code (type)
 336          //        Int_t parent = itsHit->GetParticle()->GetFirstMother();
 337           //partcode = itsHit->GetParticle()->GetPdgCode();
 338
 339    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
 340    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 341
 342            /*
 343            px = itsHit->GetPXL();  // momentum for the hit
 344            py = itsHit->GetPYL();
 345            pz = itsHit->GetPZL();
 346            pmod = 1000*sqrt(px*px+py*py+pz*pz);
 347            */
 348
 349            pmod = itsHit->GetParticle()->P(); // momentum at the vertex
 350            pmod *= 1.0e+3;
 351
 352            // cout<<"track,partcode,pmod,hitprim ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<endl;
 353
 354          if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 355                                                   // at p < 6 MeV/c
 356
 357   // find the primery particle path length in Si and pmod (MeV/c)
 358           if((hitstat == 68 || hitstat == 33) && hitprim == 1)  {
 359              pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
 360              PmodRecP = itsHit->GetParticle()->P();
 361              PmodRecP *= 1.0e+3;
 362              // cout<<"path in Si="<<pathInSi<<endl;
 363
 364          Float_t Signal = signalRP*(300/pathInSi)/38;
 365             if(Signal<0.5) {
 366               //cout<<"track,partcode,pmod,hitprim,hitstat,Signal,dEn ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<","<<hitstat<<","<<Signal<<","<<dEn<<endl;
 367             }
 368           }
 369         } // hit loop
 370
 371         if(RecPointPrim == 1) {
 372
 373           //cout<<" SDD/SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
 374
 375          Int_t xpartcode=gAlice->Particle(TrackRecPoint)->GetPdgCode();
 376
 377          if(mod < first2) {            // SDD
 378           signalRP *= (280./pathInSi); // 280 microns of Si thickness
 379           signalSDD->Fill(signalRP);
 380           signalSDDmip->Fill(signalRP/280.); // ADC/280 -> mips
 381          InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/280. ,
 382                (Float_t)PmodRecP,(Int_t)xpartcode);
 383           //cout<<" SDD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
 384           if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SDD ="<<pathInSi<<endl;
 385
 386           //if(signalRP/280 < 5) cout<<" small signalSDDmip, path ="<<signalRP/280<<","<<pathInSDD<<endl;
 387          }else{                        // SSD
 388           signalRP *= (300/pathInSi); // 300 microns of Si thickness
 389           //cout<<" SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
 390           signalSSD->Fill(signalRP);
 391           signalSSDmip->Fill(signalRP/38.); // ADC/38 -> mips
 392           InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/38. ,
 393                 (Float_t)PmodRecP,(Int_t)xpartcode);
 394           if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SSD ="<<pathInSi<<endl;
 395
 396           if(signalRP/38 < 0.5) cout<<" small signalSSD (mip), path, module ="<<signalRP/38<<","<<pathInSi<<mod<<endl;
 397          }
 398         } // primery particle
 399        } // pnt, recpoint loop
 400
 401        // dEdX hit analysis for the SDD/SSD
 402
 403          Int_t track = -100;
 404          Int_t trackPrev = -100;
 405          parent = -100;
 406          Int_t parentPrev = -100;
 407          Int_t flagprim = 0;
 408          Int_t TrackPart;
 409          epart = 0;
 410          pathInSi = 1.0e+7;
 411          Float_t PmodPart;
 412          zhit0 = 1.0e+7;
 413          xhit0 = 1.0e+7;
 414          yhit0 = 1.0e+7;
 415
 416          /*
 417          if(mod <= 259) {
 418          cout<<" SDD hits: nhits ="<<nhits<<endl;
 419          }else{
 420          cout<<" SSD hits: nhits ="<<nhits<<endl;
 421          }
 422          */
 423
 424        for (hit=0;hit<nhits;hit++) {
 425
 426          itsHit   = (AliITShit*)Mod->GetHit(hit);
 427          if(hit>0) itsHitPrev   = (AliITShit*)Mod->GetHit(hit-1);
 428          hitprim = 0;
 429          Int_t hitprimPrev = 0;
 430          track = itsHit->GetTrack();
 431          if(hit > 0) Int_t trackPrev = itsHitPrev->GetTrack();
 432          dray = 0;
 433          hitstat = itsHit->GetTrackStatus();
 434
 435           zhit = 10000*itsHit->GetZL();
 436           xhit = 10000*itsHit->GetXL();
 437           yhit = 10000*itsHit->GetYL();
 438           Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
 439
 440           //              cout<<"New hit,lay,hitstat,yhit,ehit ="<<hit<<","<<layer<<","<<hitstat<<","<<yhit<<","<<ehit<<endl;
 441           //  cout<<"befor: track,trackPrev ="<<track<<","<<trackPrev<<endl;
 442
 443            parent = itsHit->GetParticle()->GetFirstMother();
 444            if(hit>0) Int_t parentPrev = itsHitPrev->GetParticle()->GetFirstMother();
 445            partcode = itsHit->GetParticle()->GetPdgCode();
 446
 447    //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
 448    //  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 449
 450            px = itsHit->GetPXL();
 451            py = itsHit->GetPYL();
 452            pz = itsHit->GetPZL();
 453            pmod = 1000*sqrt(px*px+py*py+pz*pz);
 454
 455            //  cout<<"partcode,pmod,parent,parentPrev,hitprim,flagprim ="<<partcode<<","<<pmod<<","<<parent<<","<<parentPrev<<","<<hitprim<<","<<flagprim<<endl;
 456
 457
 458          if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 459                                                   // at p < 6 MeV/c
 460
 461
 462           if(parent < 0 && parent > -90) hitprim += 1;
 463           if(parentPrev < 0 && parentPrev >-90) hitprimPrev += 1;
 464          // hitprim=1 for the primery particles
 465           //      if(hit==0&&hitprim==1&&hitstat==66) {
 466           if(hitprim==1&&hitstat==66) {
 467             zhit0 = zhit;
 468             xhit0 = xhit;
 469             yhit0 = yhit;
 470           }
 471
 472           if((hitstat == 68 || hitstat == 33) && hitprim == 1) {
 473              pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
 474              TrackPart = track;
 475              PmodPart = itsHit->GetParticle()->P();
 476              PmodPart *= 1.0e+3;
 477              //  cout<<"xhit0,yhit0,zhit0,pathInSi ="<<xhit0<<","<<yhit0<<","<<zhit0<<","<<pathInSi<<endl;
 478           }
 479
 480           if(hitprim == 0) track = parent;
 481           if(hitprimPrev == 0) trackPrev = parentPrev;
 482           //  cout<<"after: track,trackPrev ="<<track<<","<<trackPrev<<endl;
 483
 484           if(hit > 0 && track == trackPrev) epart += ehit;
 485           if((hit > 0 && track == trackPrev) && (hitprim==1)) flagprim +=1;
 486
 487           //  cout<<"new hitprim, flagprim ="<<hitprim<<","<<flagprim<<endl;
 488
 489           if((hit > 0 && track != trackPrev) || (hit == nhits-1)) {
 490             if(flagprim > 0) {
 491               if(layer > 2 && layer < 5) {
 492                 if(epart < 40) cout<<" small SDD dedx ="<<epart<<endl;
 493                 if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSDD ="<<pathInSi<<endl;
 494                 if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) {
 495                  dedxSDD->Fill(epart);
 496                  epart *= (280/pathInSi);
 497                  dedxSDDnorm->Fill(epart);
 498                  dedxSDDmip->Fill(epart/75);
 499                  //InPID(pid,(Int_t)nev,(Int_t)track, epart/75. ,(Float_t)pmod,(Int_t)partcode);
 500                  pathSDD->Fill(pathInSi);
 501                  //  cout<<"SDD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl;
 502                 }
 503               }
 504               if(layer > 4) {
 505                 if(epart < 40) cout<<" small SSD dedx ="<<epart<<endl;
 506                 if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSSD ="<<pathInSi<<endl;
 507                 if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) {
 508                  dedxSSD->Fill(epart);
 509                  epart *= (280/pathInSi);
 510                  dedxSSDnorm->Fill(epart);
 511                  dedxSSDmip->Fill(epart/80);
 512                  //InPID(pid,(Int_t)nev,(Int_t)track, epart/80. ,(Float_t)pmod,(Int_t)partcode);
 513                  pathSSD->Fill(pathInSi);
 514                  //cout<<"SSD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl;
 515                 }
 516               }
 517               flagprim=0;
 518               epart = 0;
 519             }else{
 520               //cout<<"No!, epart for the secondary"<<endl;
 521               epart = 0;
 522             }
 523           }
 524        } // SDD/SSD hit loop
 525     } // module loop
 526     //cout<<"negtrSDD, negtrSSD ="<<negtrSDD<<","<<negtrSSD<<endl;
 527    } // event loop
 528
 529    cout<<endl;
 530
 531    TFile fhistos("dEdX_his.root","RECREATE");
 532
 533    dedxSDD->Write();
 534    dedxSDDnorm->Write();
 535    dedxSDDmip->Write();
 536    dedxSSD->Write();
 537    dedxSSDnorm->Write();
 538    dedxSSDmip->Write();
 539    signalSDD->Write();
 540    signalSDDmip->Write();
 541    signalSSD->Write();
 542    signalSSDmip->Write();
 543    pathSDD->Write();
 544    pathSSD->Write();
 545
 546    fhistos.Close();
 547    cout<<"!!! Histogramms and ntuples were written"<<endl;
 548    // ------------------------------------------------------------
 549
 550    TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
 551    c1->Divide(2,2);
 552
 553
 554 //     c1->cd(1);
 555 //     gPad->SetFillColor(33);
 556 //     signalSDD->SetFillColor(42);
 557 //     signalSDD->Draw();
 558
 559 //     c1->cd(2);
 560 //     gPad->SetFillColor(33);
 561 //     signalSDDmip->SetFillColor(42);
 562 //     signalSDDmip->Draw();
 563
 564 //     c1->cd(3);
 565 //     gPad->SetFillColor(33);
 566 //     dedxSDDnorm->SetFillColor(46);
 567 //     dedxSDDnorm->Draw();
 568
 569 //     c1->cd(4);
 570 //     gPad->SetFillColor(33);
 571 //     dedxSDDmip->SetFillColor(46);
 572 //     dedxSDDmip->Draw();
 573
 574    c1->cd(1);
 575    gPad->SetFillColor(33);
 576    signalSSD->SetFillColor(42);
 577    signalSSD->Draw();
 578
 579    c1->cd(2);
 580    gPad->SetFillColor(33);
 581    signalSSDmip->SetFillColor(42);
 582    signalSSDmip->Draw();
 583
 584    c1->cd(3);
 585    gPad->SetFillColor(33);
 586    dedxSSDnorm->SetFillColor(46);
 587    dedxSSDnorm->Draw();
 588
 589    c1->cd(4);
 590    gPad->SetFillColor(33);
 591    dedxSSDmip->SetFillColor(46);
 592    dedxSSDmip->Draw();
 593
 594
 595    /*
 596     c1->cd(1);
 597    gPad->SetFillColor(33);
 598    dedxSDD->SetFillColor(42);
 599    dedxSDD->Draw();
 600
 601    c1->cd(2);
 602    gPad->SetFillColor(33);
 603    pathSDD->SetFillColor(42);
 604    pathSDD->Draw();
 605
 606    c1->cd(3);
 607    gPad->SetFillColor(33);
 608    dedxSSD->SetFillColor(46);
 609    dedxSSD->Draw();
 610
 611    c1->cd(4);
 612    gPad->SetFillColor(33);
 613    pathSSD->SetFillColor(46);
 614    pathSSD->Draw();
 615    */
 616
 617 //   if(!pid)pid->Tab();
 618    cout<<"END  test for clusters and hits "<<endl;
 619
 620 }
 621 //============================ InPID() =========================================
 622 int totpid;
 623 void InPID (AliITSPid *pid=0,int     nev,
 624                              Int_t   track,
 625                              Float_t signal,
 626                              Float_t pmom,
 627                              Int_t   pcode)
 628 {
 629 // Includes recp in PID  table.
 630     if(pid){
 631       //cout<<" In pid: track,sig,pmod,pcode="<<track<<","<<signal<<","<<pmom<<","<<pcode<<endl;
 632     if( abs(pcode)==321 || abs(pcode)==211 ||abs(pcode)==2212 ||abs(pcode)==11  )
 633         {
 634         pid->SetEdep(10000*nev+track,signal);
 635         pid->SetPmom(10000*nev+track,pmom);
 636         pid->SetPcod(10000*nev+track,abs(pcode));
 637         }
 638     }
 639     totpid++;
 640 }
 641 //=======================
 642
 643
 644
 645
 646
 647